Finite-size scaling of directed percolation in the steady state 
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Recently, considerable progress has been made in understanding finite-size scaling in equilibrium 
systems. Here, we study finite-size scaling in non-equilibrium systems at the instance of directed 
percolation (DP), which has become the paradigm of non-equilibrium phase transitions into absorb- 
ing states, above, at and below the upper critical dimension. We investigate the finite-size scaling 
behavior of DP analytically and numerically by considering its steady state generated by a homoge- 
neous constant external source on a d-dimensional hypercube of finite edge length L with periodic 
boundary conditions near the bulk critical point. In particular, we study the order parameter and its 
higher moments using renormalized field theory. We derive finite-size scaling forms of the moments 
in a one-loop calculation. Moreover, we introduce and calculate a ratio of the order parameter 
moments that plays a similar role in the analysis of finite size scaling in absorbing nonequilibrium 
processes as the famous Binder cumulant in equilibrium systems and that, in particular, provides a 
new signature of the DP universality class. To complement our analytical work, we perform Monte 
Carlo simulations which confirm our analytical results. 

PACS numbers: 64.60.Ak, 05.70.Jk, 64.60.Ht 



I. INTRODUCTION 

Critical phenomena like second order phase transitions 
are characterized by singularities of various quantities at 
the transition point (e.g. the specific heat, susceptibil- 
ity, correlation length) . These singularities are described 
by power-laws governed by critical exponents. Studying 
the phase transition of a given system, one usually tries 
to identify the set of critical exponents which in con- 
junction with certain universal scaling functions charac- 
terizes the present universality class. Powerful analyti- 
cal and numerical techniques have been developed to ac- 
complish this task. Analytical investigations of universal 
quantities allow to address infinite system sizes but they 
are usually feasible only if one uses involved approxima- 
tions such as the diagrammatic perturbation expansions 
of renormalized field theory. Using numerical techniques 
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like Monte Carlo simulations or transfer matrices calcu- 
lations one can avoid such approximations, however, the 
data is limited to finite systems sizes. Therefore, finite- 
size scaling (FSS) is widely used to extrapolate to the 
behavior of infinite systems. In particular, FSS is an 
efficient method to determine critical exponents and cer- 
tain universal scaling functions, and therefore, it often 
allows to identify the universality class (see Refs. [I, 2] 
for reviews). According to the phenomenological FSS 
theory [3] , finite system sizes L result in a rounding and 
shifting of the critical singularities. It is assumed that 
finite-size effects in isotropic systems are controlled suffi- 
ciently close to the critical point by the ratio i/^oo, where 
^oo is the spatial correlation length of the infinite system. 
Approaching the transition point, this correlation length 
diverges as ^oo r~'^ , where r oc |t — rd measures the 
deviation of a temperature-like control parameter t from 
its critical point value Tc, and where is the critical expo- 
nent of ^oo • Finite-size effects decrease with increasing L 
and are negligible for L ^ ^oo, i-e., for L^^^r ^ 1, in sys- 
tems with periodic boundary conditions, true short range 
interactions, and without Goldstone modes. Otherwise, 
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they are relevant, i.e., rounding and shifting effects occur 
when L < ^^a. It is well known that in equilibrium the 
hypothesis of the fundamental role of the ratio L/S,oo is 
valid only below the so-called upper critical dimension dc 
(see [4] for a recent review). Above dc, mean field theo- 
ries provide exact results for the critical exponents and 
the scaling functions. However, usual FSS fails above dc 
because certain parameters, which are irrelevant in the 
sense of the renormalization group, become dangerously 
irrelevant for d > dc [5] . Dangerous irrelevant parameters 
affect the scaling behavior qualitatively and furthermore 
cause the breakdown of hyperscaling laws which connect 
the critical exponents to the spatial dimensions d. Inves- 
tigations of this breakdown of usual finite scaling date 
back to the work of Brezin and Zinn- Justin [6-8] . For the 
case of periodic boundary conditions, Brezin and Zinn- 
Justin introduced an analytic technique which makes it 
possible to perform calculations of size dependent uni- 
versal scaling functions. This method exploits the fact 
that the so-called lowest or zero mode is distinguished 
in the sense that in perturbation theory it becomes crit- 
ical before the higher modes do and that, therefore, the 
latter modes can be traced out perturbatively and fully 
neglected above dc- 

It must be emphasized that a meaningful, quasi- 
universal analytical study of finite-size effects is possible 
only in a regime where 1 /L <C 1 and r <C 1 , where it is 
understood that L and r are measured in terms of suit- 
able non-universal amplitudes. Outside this regime, in 
particular, if L becomes smaller, finite-size effects will be 
blurred by the effects of variables that are irrelevant with 
respect to the corresponding bulk universality class. For 
L = 0(1), analytic approaches are ultimately hopeless. 
Above dc, the strongest irrelevant effects stem from the 
usual coupling constant (in the following denoted g) of 
the non-harmonic term in the field theoretic functional, 
which is relevant below dc, but which is dangerously ir- 
relevant above dc- 

After controversial discussions of the zero-mode theory 
and the influence of the higher modes (see e.g. [4, 9, 10] 
and references therein) the problem was recently resolved 
by Chen and Dohm [10], and convincing agreement be- 
tween numerical data and field theoretical results was 
achieved [11]. Chen and Dohm showed that even above 
the upper critical dimension dc the higher modes play 
an essential role. The following three points summarize 
key findings: (i) The higher modes induce a shift of the 
critical value of the control parameter proportional to 
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Figure 1: Scaling regions (schematically) above dc where cor- 
rections to the lowest mode approximation resulting from 
higher modes are essential (I and III) or negligible (II). 

g2j2-d^ where g is the dangerously irrelevant coupling 
constant, cf. region I in Fig. 1. This shift is crucial for 
the correct interpretation of simulations, (ii) The in- 
fluence of the higher modes is essential for the correct 
description of the exponential decrease of the finite size 
effects approaching the infinite volume limit, cf. region 
III in Fig. 1. Points (i) and (ii) suggest that the correc- 
tions induced by the higher modes can be neglected only 
in the region 

g2^2^d/2 ^ (1.1) 

cf. region II in Fig. 1 [12]. (iii) Chen and Dohm shed light 
on the fact that analytical methods using a hard momen- 
tum cutoff, which is well known to be equivalent to long 
range interactions, induce a wrong algebraic decrease of 
finite size effects. Hence, the widely used Fisher- Wilson 
momentum shell, like any other hard-cutoff renormaliza- 
tion procedure, is incompatible with the exponentially 
decreasing crossover to the infinite volume limit. 

Compared to the equilibrium situation, much less is 
known in the case of non-equilibrium phase transitions. 
This motivates us to discuss in the paper at hand FSS 
in non-equlibrium phase transitions at the instance of di- 
rected percolation (DP). Due to its robustness and ubiq- 
uity (including critical phenomena in physics, biology, 
dynamics of populations, epidemiology, as well as au- 
tocatalytic chemical reactions) DP is recognized as the 
paradigm of non-equilibrium phase transitions into ab- 
sorbing states (see [13-15] for a recent review on ab- 
sorbing state transitions, and [16] for a recent review 
on renormalized field theory applied to percolation pro- 
cesses) and, although an exact analytical solution is still 
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lacking, DP plays a role for non-equilibrium phase tran- 
sitions comparable to that of the Ising model for equilib- 
ria. Previous studies of finite-size scaling of percolation 
processes by one of us and coworkers [17] focused on the 
absorbing phase below dc- Here our scope is different: we 
are interested in finite-size properties of the steady state 
below, at and above dc = 4. To be specific, we study 
for these dimensions the scaling behavior of finite DP 
systems in the active phase which is maintained by a ho- 
mogeneous external source. Using Reggeon field theory 
(RFT) [19], the generic field theoretic description of the 
DP universality class [20-24], we derive finite-size scal- 
ing exponents and universal scaling functions for periodic 
boundary conditions. For d > 4, we demonstrate that 
the usual phenomenological FSS theory for DP has to 
be modified, analogous to what we have discussed above 
for the equilibrium case, in order to describe the scaling 
behavior within the mean field regime. We show that 
the correct scaling variable in the strong finite size re- 
gion L <^ ^oo oc is proportional to L'^^^r, and that 
corrections, which are controlled by an expansion in a 
variable v oc gL^^"^^^, become essential only if this vari- 
able goes to zero. 

Compared to the equilibrium case, an additional con- 
ceptual problem arises in dynamics: to obtain analytical 
results for the finite-size scaling functions, one is forced 
to perform a Markovian approximation of the dynam- 
ics of the lowest mode. Therefore, our analytical results 
are restricted to the strong finite-size region. Outside 
this region the Markovian approximation leads to a de- 
scription of the crossover to the infinite-volume limit by 
algebraically decreasing correction terms instead of the 
correct exponentially decreasing ones, even if we include 
the one-loop corrections arising from the higher modes. 
We explicitly demonstrate this failure of the Markovian 
approximation via a perturbation calculation of the cor- 
relation function. 

In the strong finite-size region, we observe, when in 
the region near the critical point the shift induced by 
the higher modes is taken into account, convincing quan- 
titative agreement between the lowest mode finite-size 
analysis and our numerical results. For d < 4 and d = A 
we calculate e-expansions and logarithmic corrections, re- 
spectively, for various quantities, focussing, in particular, 
on an universal ratio of order parameter-moments. For 
d < 4, in addition, we perform simulations which clearly 
underscore that this ratio is a universal signature of the 
DP class. Brief account of parts of the work presented 



here has been given previously in Ref. [18]. 

The outline of our paper is as follows. In Sec. II we 
briefiy review RFT as the field theoretic model of choice 
for the DP universality class. We derive the effective 
response functional, i.e., the dynamic free energy of the 
homogeneous (lowest) mode. In Sec. Ill we calculate this 
dynamic free energy in a 1-loop Markovian approxima- 
tion. In section IV we derive finite-size scaling forms for 
spatial dimensions above the upper critical dimension. 
The steady state solution of the Fokker-Planck equation 
which correspondents to the effective response functional 
yields all moments of the homogeneous mode in scaling 
form. In Sec. V, we compare our analytical results with 
numerical results stemming from our Monte Carlo simu- 
lations. In Sec. VI we study the crossover to mean field 
theory in the infinite volume limit. In Sec. VII we con- 
sider finite size effects in the steady state for spatial di- 
mensions below the upper critical dimension. We apply 
the renormalization procedure to our 1-loop results, and 
we derive universal values of the aforementioned ratio of 
order-parameter moments in an e-expansion. The ana- 
lytic estimates which follow from this expansion are com- 
pared with the numerical results. In Sec. VIII, we study 
finite size effects right at the upper critical dimension. 
We calculate logarithmic corrections to various quanti- 
ties including our momenta ratio. Concluding remarks 
are given in Sec. IX. An appendix contains a brief pre- 
sentation of the properties of some functions fundamen- 
tal to finite-size scaling in DP. For the convenience of the 
reader, we will provide at the beginning of the main sec- 
tions short summaries of their respective contents and we 
point out to their most important formulas. 

II. REGGEON FIELD THEORY AND THE 
EFFECTIVE RESPONSE FUNCTIONAL 

We start our analysis by deriving an effective response 
functional for the zero mode. One of the main findings 
of this section is that the distance r from the bulk crit- 
ical point in this effective theory is given by Eq. (2.7). 
Equation (2.11) summarizes our result for the effective 
response functional. 

It has been known for a long time that the DP uni- 
versality class is well represented by RFT. For a recent 
overview on the field theories of percolation processes and 
the derivation of the underlying minimal models from 
basic principles see [16]. RFT, based originally upon a 
non-hermitean Hamilton-operator [19], is equivalent to a 
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Langevin-description of a minimal DP process, the so- 
called Gribov process [20]. After reduction to the rel- 
evant terms, the stochastic equation of motion of this 
DP process may be written in the form of the Langevin 
equation (in the Ito interpretation) [23, 24]: 

X-^dts{r,t) ^ - V2+r+|s(r,t) s(r, t) + + ((r, t). 

(2.1) 

Here, the activity field s(r, t) > 0, which is proportional 
to the density of active particles (agents) on a mesoscopic 
(coarse grained) scale, is the order parameter field of the 
non-equilibrium phase transition. The diffusional term 
represents the isotropic spreading of activity. The con- 
trol parameter of the transition is r, and denotes its 
critical value. In the infinite volume limit, a finite pos- 
itive particle density occurs below the transition point 
(r < Tc) whereas the absorbing vacuum state (s ~ 0) is 
approached above the transition point if the source term 
h > (which can be implemented in simulations, e.g., as 
a spontaneous particle creation process [25]) is absent. 
In a finite system, the absorbing state is inevitably ap- 
proached even for t < Tc, if h = 0. However, it can be 
shown [27] that the logarithm of the relaxation time to 
the absorbing state increases proportional to the system 
volume in the active phase below Tc. A and g denote the 
kinetic and coupling constants, respectively. C: finally, 
represents the noise which accounts for fluctuations of 
the particle density. All universal properties of the DP 
universality class arc captured by the minimal model, 
Eq. (2.1), provided the noise Ci^i^) is a Gaussian ran- 
dom variable with zero mean and correlator given by [26] 



considered stochastic process. The dynamic response 
functional of the Gribov process (2.1) is given by [23, 24] 



C(r, t)C(r', t') = A-\g' s(r, t)5(r - r')6{t - t'). (2.2) 

Note, that only an absorbing noise with a correlator that 
comprises at least a term linear in the field s ensures 
that the systems is trapped in the absorbing state with 
a continuously decreasing survival probability [27]. A 
form of the noise proportional to (multiplicative noise) 
results in a survival probability which is strictly 1 for all 
finite times. 

Renormalization group techniques have been applied 
quite successfully to determine the critical exponents and 
universal scaling functions of DP [16, 21-24, 27-29]. In 
the framework of field theory, a path integral formulation 
of stochastic processes is more useful than their Langevin 
equations. In the path integral formulation, correlation 
functions and response functions can be determined by 
calculating path integrals with weight exp {—J) [30, 31], 
where the dynamic response functional describes the 



J[s,s]= / d'^rdu{s[A-iat + (T-V2) + |(s-S) 



(2.3) 

where s(r,t) denotes the purely imaginary response field 
conjugated to the Langevin noise field. The functional 
is invariant under time reversal (in RFT usually called 
rapidity reversal) 



~s{Y,t) 



- s(r, -t) 



(2.4) 



as long as the (symmetry breaking) field h vanishes. This 
symmetry is spontaneously broken in the active phase 
below the transition point. In general, the time reversal 
invariance of the minimal model is merely an asymptotic 
symmetry of systems belonging to the DP universality 
class. Note, however, that this symmetry is exact for 
bond DP. 

It is worth noting, that the original RFT [19] is based 
on a bosonic annihilation-creation formalism in which s is 
related to the annihilation operator, and s to the creation 
operator (for a recent review over this master-equation 
approach see [32]). Hence, as described in [16], the orig- 
inal RFT and the fiuctuating field theory based on the 
Gribov process (2.1), where s is proportional to a real 
positive density, are formally different. Note that the 
bosonic theory leads to an additional noise term propor- 
tional to {ssY in the functional J , Eq. (2.3), with pos- 
itive sign. Such noise terms result from anticorrclating 
more-particle annihilation reactions, and are typical for 
diffusion-limited reactions. However, for DP, which is not 
a diffusion-limited reaction system, this noise term is ir- 
relevant. Hence, both formalism, the Langevin- and the 
master-equation approach, produce the same perturba- 
tion series which leads via the renormalization group to 
the same universal asymptotic behavior. Note, however, 
that in the bosonic formalism s and s arc constructed as 
complex fields with ss real and positiv. Thus, after delet- 
ing the irrelevant noise term, the functional integration 
with the weight exp (— i/) is a priori mathematically ill- 
defined. Ciafaloni and Onofri [33] have shown more than 
25 years ago that in this case the only correct support for 
integration over s and s is, respectively, the real positive 
axis and the full imaginary axis (see also the appendix 
of [24] for a corresponding quasi-canonical transforma- 
tion of the fields). The upshot is that only the Langevin 
equation formalism offers a mathematically correct in- 
terpretation of the functional integral. Because we must 
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use at least parts of the weight exp {—J) as it stands, 
without applying perturbation theory, as a probabihty 
measure, this interpretation is of greatest importance. 

Using standard techniques known from equihbrium [6- 
8], one can extract from an effective response func- 
tional for the zero mode, which then can be used to cal- 
culate size-dependent universal scaling functions as well 
as the involved critical exponents. To follow this route, 
let us consider DP in a finite cubic geometry of linear 
size L with periodic boundary conditions and expand s 
and s in plane waves. 



s{v,t) 



q 



e''i--s(q,i). 



(2.5) 



and likewise for s. Each component of the wavevector q 
takes on discrete values, viz. multiples of 2t: / L including 
zero. When dealing with summations over q, one has to 
bear in mind that path-integrals based on the response 
functional (2.3) are well-defined only if an appropriate 
regularization of the diverging UV behavior is applied. 
In principle, there are different options for choosing a 
regularization procedure. As discussed in the introduc- 
tion and as can be easily checked by applying the Euler- 
McLaurin summation formula, a hard momentum cutoff 
(support of the modes only for momenta with |q| < A) 
is inappropriate for studying FSS, since a hard cutoff in- 
duces non-physical long-range correlations in real space 
[8] which contaminate the finite size calculations [10]. 
Lattice regularization, where the system is placed on a 
discrete lattice instead of spatial continuum, is the most 
physical one. Moreover, this regularization is in closest 
contact to simulations. However, lattice-regularization 
replaces the Laplacian by the lattice difference-operator. 
Thus, analytical calculations become very complicated. 
In the following, we will use (implicitly) a soft cutoff 
procedure, i.e., we will include a factor exp(— q^A^) in 
all summations over wave vectors q, followed by dimen- 
sional regularization and the limit K ^ oo. One can 
show that this procedure is equivalent in the scaling re- 
gion to lattice regularization as long as one concentrates 
on universal quantities. At this point, a word of caution 
is in order. If very small lattices are considered, it may 
be more appropriate to use lattice regularization [10]. In 
the following, we will ignore very small lattices in our 
analytic considerations because for these lattices one has 
to expect many non-universal corrections. 

The Fourier transformed propagator of the perturba- 
tion theory about the saddle-point of the path-integrals 
(mean-field theory) based on the response functional 



Eq. (2.3) is given by 

Go(q,a;) 

with 



iuj/\ + 7' + q2 



r = (t ~Tc) + M 



(2.6) 



(2.7) 



measuring the distance to the critical point. Here M = 
g{s), with the expectation value (s) determined by the 
condition that tadpoles are excluded in the diagram- 
matical perturbation expansion. In mean-field theory, 
r = \J t'^ + 2gh and Tc = 0. Hence, for small frequen- 
cies u! in the finite-size limit, w = r(L/27r)^ ^ 1, the 
zero mode with q = separates from the higher modes 
and leads to infrared divergencies in perturbation theory. 
Therefore, functional integrals of the zero mode must be 
calculated exactly, and cannot be handled by perturba- 
tion theory [6-8]. Perturbation theory can be used, how- 
ever, as a tool for the functional integration of the higher 
modes. As we will discuss in detail later on, the Gaus- 
sian fluctuations of the higher modes have a significant 
influence on the scaling functions describing the crossover 
from If w 1 to w 3> 1, as well as the behavior near the 
bulk critical point w « 0. Nonetheless, mean-field theory 
should be correct for d > 4 in the bulk limit w ^ oo. 

Following [17], we construct an effective response func- 
tional for the zero-mode by separating the homogenous 
mode ^{t) from its orthogonal complements \E'(r,t) via 
setting 



gs(r,t) = $(<) + «'(r,<) 



(2.8) 



with <I>(t) ™ gL^"^ J d'^rs{r,t) and likewise for s. This 
leads to a decomposition of the action, J — J0 + J1+J2, 
with 



Jo = \g-^L'' j dt{l>[A-iat + T + i($-$) 



(2.9) 



where H = gh^ and 



1 



($ - $) 



^1 



')}■ 



Ji^X^ / dVd< «'(«'- *)*. 



(2.10a) 
(2.10b) 



We have included the coupling constant g in the defi- 
nition of the fields <!>, $, , ^ to disentangle the two 
different roles of which on the one hand serves as the 
loop-order generating parameter of the perturbation the- 
ory around the mean-field (Landau) approximation, and 
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on the other hand is a scale factor of the fields. This last 
role is what makes g a "dangerous" irrelevant variable, 
as alluded to in the introduction. Finally, we eliminate 
^ and via functional integration. This leads to 

E[l>,$] = Jo -In J V[^, *] exp(- Ji - J2) (2.11) 

as our effective response functional for the homogeneous 
mode. In the following, we will also refer to E as our 
dynamic free energy. 

The zero-loop approximation E sa j7o is known as the 
lowest mode approximation of finite-size scaling [7, 8] . As 
we move along, we will show that, for > 4, this lowest 
mode theory is modified outside the lowest-mode region 
^ ^ i2-d/2^ pjg_ (j)^ one-loop 

(Gaussian) contributions arising from the higher modes. 



III. DYNAMIC FREE ENERGY IN THE 
ONE-LOOP EXPANSION 

In this section, we calculate the dynamic free energy 
E to 1-loop order in perturbation theory. Key formulas 
of this section are Eqs. (3.11) and (3.15) which give the 
finite-size scaled version of, respectively, the distance r 
from the bulk critical point and the dangerously irrele- 
vant coupling constant g. Central to our discussions to 
follow is the 1-loop dynamical free energy (3.5) in con- 
junction with our 1-loop results for the parameters ap- 
pearing in it, Eq. (3.14). 

To 1-loop order. J2 does not contribute and hence can 
be neglected. jTi contributes via the propagator 



G(t,t';q) = 0{t - t') cxp 



A(r + q2)(<-<') 



+ \ dt" {^{t") - ^{t")) 



(3.1) 



of the higher modes, which is determined by the bilinear 
part in the fields '5, ^' of J7i. Gaussian integration yields 
readily 

-In /l?[«',*]exp(-Ji) = yJI j j dtdt' ^{t) 

X G(t,t';q)2$(t') + C'(($$)2). (3.2) 

For the time being, let us concentrate on the region w <C 
1. Then, the typical time-dependence of the zero-mode 
shows slowing down in comparison to the higher modes 
leading to Markovian behavior of the zero-mode. Thus, 



we can approximate $(i') in Eq. (3.2) by <^{t)-{t-t')^{t) 
and the propagator simplifies to 



G(<,t';q) = e{t-t') exp 



-A( 



(3.3) 

Note that this Markovian approximation does not any 
longer allow a correct description of the crossover from 
the finite size to the infinite volume behavior. If one in- 
correctly takes w ^ 1 in the results following from this 
approximation one gets algebraically decreasing correc- 
tion terms describing the crossover to the infinite volume 
limit. This crossover is qualitatively wrong because the 
corrections must be exponentially decreasing. We will 
discuss this shortcoming of the Markovian approximation 
in Sec. VI, where we calculate the steady state correlation 
function for w ^ 1 in a 1-loop calculation. 

After application of the Markovian approximation the 
residual time integration of t' can be done. We obtain 



In j exp(-^i) 



4(t + $ - $ + q2 



8 r 



$ - $ + q2) 



(3.4) 



retaining only terms of the form already appearing in 
Joi i-e., neglecting fourth-order terms in $ and $. These 
higher order monomials lead to corrections of higher or- 
der in L^^ as the retained ones [7, 8]. Subsequently, we 
expand the denominators in Eq. (3.4) in <!> and $ about 
their mean values ($) = M and (l>) = 0. Note that this 
procedure provides strictly positive denominators even in 
the case Tc — r > {2tt/LY and also that we can include 
the bulk critical value of the control parameter r in 
the denominators of Eq. (3.4) since is of order g^ . Re- 
calling definition (2.7), we finally obtain from Eqs. (2.9) 
(2.11) and (3.4) that 

/ 



E[$,$] 



dt <^ $ 



The parameters k, f and / are given by 

,2 



k{r) 
T{r) 
fir) 



1 



l-^52(r) 



4 



1 



\S2ir)^ 



with 5*; defined by 

Si{t) = L- 



q#0 



(3.5) 

(3.6a) 
(3.6b) 
(3.6c) 

(3.7) 
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As mentioned earlier, all sums over wavevectors must be 
regularized appropriately; it is understood that in actual 
calculations these sums are augmented by a soft cutoff 
factor cxp(— q^A^). In the infinite-size limit, L oo, 
the sums Si{r) tend to the integrals 



sr{r) 



1 



(3. 



q Ir + 



with /q • ■ • = (Stt)"'' J d'^q.... The bulk critical point is 
then given in 1-loop approximation by 



■9 Qoo,n^ y f 1 



(3.9) 



For the steps to follow, it is useful to introduce the 
differences 



Ai{r)^Sr(r)-Siir) 
with the scaling variable 



L 



2l-d 



D^'\w) (3.10) 



w 



27r 



(3.11) 



The D(') are functions of this scaling variable given by 

rf/2 



Jo 



(3.12) 



where A{t) = 1 + 2 ^^^^ exp(-n2i) = (7r/0^/=^A(7rVt). 
Some important properties of the functions D^^^w) are 
discussed in the appendix. With help of the differ- 
ences (3.10), we can express the parameters appearing 
in the dynamic free energy after some rearrangements 
as, 



fc(r) 
f(r) 



[14 


f ' 1 






/ ' 1 


iq (r + q2)2j 



+ ^A2(r), (3.13a) 
(r-T-e) 
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4 Jq<fir + (f)^ 

^ [A2(r) (r - Te) - Ai(r) - rA2(r)] , (3.13b) 



1-^ 



1 



2 L (r + q2)2 



^A2(r), (3.13c) 



where we have neglected terms of order g* . The integrals 
over wavevectors in Eqs. (3.13) lead to IR singularities 
for spatial dimensions d < 4 if r ^ 0. These singularities 
must be treated by the renormalization group [17]. We 



will return to the cases d < 4 and d = 4 in Sees. VII and 
VIII, respectively. 

For d > 4, the integrals lead to cut-off depen- 
dent nonunivcrsal constants up to corrections of order 
j-id—i)/2^ We neglect these corrections, and include the 
nonunivcrsal constants in a rescaling of the fields 
and of the parameters T,g,H. We redefine r — Tc — )■ t 
and thus, henceforth, r = at the bulk critical point. 
Finally, we obtain for d > 4 



w] 



(2)(u;) 
f= 1 + 2v^D'^^\w) 



wD 



(^) 



(3.14a) 

(3.14b) 
(3.14c) 



where we have defined a second scaling variable 

v = A.L^-''/^. (3.15) 

Now, after having identified w and v as fundamental 
scaling variables, it is worthwhile to briefly reconsider 
the condition for the approximations that we made in 
this section. To justify the neglect of higher loop-orders 
of the perturbation expansion as well as the influence of 
other irrelevant couplings in the response functional (2.3) 
we have to assume <C 1, that means that L is suffi- 
cient large but finite. Moreover, for the application of the 
Markovian approximation, we have to assume w <C 0. 



IV. FOKKER-PLANCK EQUATION AND 
SCALING OF THE STEADY STATE 
OBSERVABLES ABOVE 

In this section we analytically derive scaling forms and 
scaling functions for steady state observables above the 
upper critical dimension. First, we identify further fun- 
damental scaling variables, namely the finite-size scaled 
control parameters given in Eq. (4.1). Moreover, we in- 
troduce finite-size scaled fields, Eq.(4.3a), and a finite 
size-scaled time, Eq. (4.3b). This leads to a finite-size 
scaled dynamic free energy, Eq. (4.5), with parameters a 
and b given in Eq. (4.7), which will play a central role as 
we move along. Then, we discuss how we can calculate 
the moments of the homogeneous density, i.e., averages 
of powers of $, with the help of a Fokker-Planck equa- 
tion, Eq. (4.11), and its stationary solution, Eq. (4.12). 
Our results for the moments of the homogeneous density 
are presented in Eq. (4.18). 
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As we move along, we shall see that, besides w and v, 
three further combinations of the four physical param- 
eters T, h, M, L and the dangerous irrelevant coupling 
constant g emerge as natural variables of the finite-size 
scaling forms of the moments of the homogeneous den- 
sity, namely 



9 



y=-L^H, 



(4.1) 



The variables w and v are related to x and z by 

w ^ v{x + z) . (4.2) 

Moreover, we introduce scaled fields ip, ip, and a scaled 
time s, 

$(t) = aip{s) , $(i) = aip{s) , (4.3a) 

Xt = /3s , (4.3b) 

with scale factors a, P given by 

a = fc-i/2.9L-'^/2, (4.4a) 

f3 = 2f-^p/^ g-^L'^^\ (4.4b) 

Incorporating these rescalings into the dynamic free en- 
ergy (3.5), we obtain 

T,[i^,ip] ^ J ds ^(p ds + a + {ip ~ (fi) ip - 6(^1 , (4.5) 

with new parameters 



b = 2kg-^hL'\ 



(4.6a) 
(4.6b) 



where g = gf. When expressed in terms the scaling 
variables, these new parameters read 



l_^i^(2)(^)+0(«4. 



1 



D^^\w)+wD^^^{w)+0{v^) 



-D^^\w)+0{v^) 



y- 



(4.7a) 
(4.7b) 



Here we have indicated 2-loop contributions and higher 
order ones, resulting from neglected irrelevant couplings 
in the response functional (2.3), by the Landau order 
symbol. The rescaled form (4.5) makes transparent an 



essential feature of the dynamic free energy S[(^,</3] : it de- 
pends only on the two parameters a and b. As a straight- 
forward consequence, all correlation and response func- 
tions of the homogeneous density (which is proportional 
to if) are universal functions of a, b, and the scaled time 
s: 

I N N+N \ 

i^^,^({s J, a, 6) = / \{^{s,) n p>{s,) \ . (4.8) 

The strict lowest-mode approximation neglects all the 1- 
loop corrections of the higher modes, that is it sets v = 0. 
Hence, we have a{v = 0) = x and b{v ~ Q) ~ y in this 
approximation. Recalling the definitions (4.3), we find 
that the correlation and response functions 

IN N+N \ f™'") 

G^,,{{U},T,h,L,g,\) = IWHU) n ^fe)) 

(4.9) 

of the homogeneous density have the finite-size scaling 
form in the zero-loop approximation 

G^,^({i«}, r, K L, g, A) - {g'/ L'')(^+^^/^ 
X Fj,j,{{L-^/^gXU/2}, 2i'^/ Vff, '^L'h/g) (4.10) 

with the universal scaling functions Fjy jy. 

Next, we determine the scaling functions 
^jv,o({0}, a, 6) including their 1-loop corrections. 
Path integrals with weight exp{—'S[(f,(p]) based on the 
dynamic free energy S[i^,iy9], Eq. (4.5), are equivalent to 
mean values taken with a probability P(ip, s\ipo)dip to 
find the process in the interval [(^s, if + dip] at time s if 
the process is started at time with ipQ. The probability 
density P{ip, s\ipo) is determined by the Fokker-Planck 
equation [17, 27] 

dsP{ip, s\lpo) = d^{[{a + - b]P{ip, s|(po)} 

+ 52{^F(^,,s|^o)}, (4.11) 

with initial condition P(iy9, 0|iy9o) = ^{'f^'^o)- In the clas- 
sification scheme of Feller [34, 35] , the infinite point ip = 
oo is a natural boundary, and therefore P(cx), s|(y5o) = 0. 
The boundary at = is a so-called exit boundary, rep- 
resenting the absorbing state as a growing J-function, if 
6 = 0. In the case 5 > 0, this boundary is regular (en- 
trance) if < 6 < 1, and natural for 5 > 1. In both cases 
it is easy to find the steady state distribution: 



2 

Po{ip) = C ip'-^ exp(-a^-^) 



(4.12) 
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where C is determined by the normalization condition 
dip Po{f) = 1. Note that in the hmit 6 — > the 
normahzation constant C goes to zero as a consequence 
of the absorbing state. In this case the only normalizablc 
stationary probability density is limt^g Poiv) = ^{'p). 

Now we fix the scaling variable z. Because a and b are 
the only parameters that our state distribution depends 
on, z enters the 1-loop correction terms, but it does not 
appear at zero-loop order. Hence, we can use here the 
strict lowest mode approximation z = 2((^)o, where the 
mean value is calculated with the steady state distribu- 
tion (4.12) with a and b set equal to the zero-loop forms 
a{v ^ 0) = X and b{v = 0) = y. This leads us to 



where iJi is a member of the set of functions defined by 
(V?') ^i{a,b), i.e.. 



Ma,b) 



dip Lp'^~^ exp(^—aip — ip^ 
rib + l)D^b^i{a) 
T{b)D^b{a) ■ 



(4.14) 



z = 2Mx,y): 



(4.13) 



Here, Da{z) are the well known parabolic cylinder (We- 
ber) functions (see, e.g., [36]) . Using the relations of 
these functions to the confluent hypergeometric (Kum- 
mer) functions M{a,(j]z) with M(a, /3;0) = 1, we have 



\/2r((l + 6)/2)A/((l 4- 6)/2, 1/2; - 2ar(l + V2)M(1 + 6/2, 3/2; aV2) 



r(6/2)M(6/2,l/2;a2/2)- 



V2ar((l 
I 



6)/2)A/((l-f &)/2,3/2;aV2) 



(4.15) 



and the recursion relation 

Ma.b) = (b + I - 2) ^i^2ia,b) - a^i.iia,b) (4.16) 

with i?o(aj b) — 1. We note two other special relations for 
later use: 



- o'/2 



m + b)/2) 
T{b/2) 



(4.17a) 



?9i (a, b) = b v/V2 e"'/^ erfc(a/V2) + 0(6^) . (4.17b) 

Here, erfc(x') denotes the complementary error function. 

Now we are finally in the position to write down a 
scaling form for the moments of the homogeneous density 
with known scaling functions. Collecting, we obtain 

(<I)^) = {kL'/gY''^^^Nia,b), (4.18) 

with universal scaling functions i^jv,o({0}, a, 6) = 
'!?7v(a,6) given to one-loop order by Eq. (4.15), or re- 
spectively, immediately following from Eq. (4.15) via the 
recursion relation (4.16). 



V. SIMULATION RESULTS ABOVE dc 

To complement our analytical calculations for d > dc, 
we have performed Monte Carlo simulations at d = 5 of 
two critical models belonging to the DP universality class 



(see [13] and references therein). Naturally, the observ- 
ables that we found best suited for our numerical work 
were not necessarily identical to those that are most con- 
venient for doing field theory. In the following, we iden- 
tify observables (ratios) that are convenient for numerical 
work, Eq. (5.4). Then, we connect these observables with 
our field theoretic results which provides us with scaling 
fimctions for these observables, Eq. (5.5). We introduce 
a ratio U, Eq. (5.6), which in a certain sense takes on 
the role in critical dynamics that the famous Binder cu- 
mulant plays in equilibrium critical phenomena. Equa- 
tion (5.7) gives our general analytical result for U. We 
derive the scaling form of the ubiquitous parameter a 
at the critical point, Eq. (5.11). This finally leads us 
to Eq. (5.12) for U, which is particularly well suited for 
comparison between theory and simulation. 

We have simulated the contact process (CP) on sim- 
ple cubic lattices of size L = 4, 8, 16 at the criti- 
cal value of the respective control parameter A, A = 
Ac = 1.13846(11), as well as the site-directed percolation 
process (sDP) implemented via a generalized Domany- 
Kinzel automaton [38, 39] on bcc lattices of linear size 
L = 8, 16, 32 at the critical value of the occupation prob- 
ability p,p = pc = 0.0359725(2) [40]. In contrast to con- 
ventional equilibrium simulation techniques, steady state 
finite-size quantities arc inaccessible for absorbing phase 



10 



transitions at zero field because, close to the transition 
point, the systems will be soon trapped in the absorb- 
ing state without chance of escape. To circumvent these 
difficulties, we perform simulations in non-zero source at 
criticality, as recently advocated in [37]. In remainder 
of this section we will present the results of our simula- 
tions and compare them to the analytic results derived 
in Sec. IV. 

Using first the lowest mode approximation without the 
1-loop corrections of the higher modes we are interested 
in the moments of the order parameter, the homogeneous 
density ^, about the absorbing state $ = 0. According 
to Eq. (4.10), we have the scaling equations 

(AiL)^^/2(ci>^) = MN{ArTiALLy/^Ahh{ALLf) . 

(5.1) 

Deviating from the conventions used in Sec. IV, we have 
here explicitly pulled the non- universal amplitudes A^, 
Ah, and A^ out of the parameters L, h and r, re- 
spectively. In accord with our analytical result (4.17a), 
we choose the normalizations Afi(0,l) = \/2/t: and 
M4(0, 1) = 3(71/2(0, 1))^. With these normalizations, we 
get for the universal finite size scaling functions defined 
in Eq. (5.1): 



Miv(0,y) 



riy/2) 



(5.2) 



where bulk criticality, r = 0, is assumed. For the order 
parameter, in particular, this leads to the modified FSS 
scaling form 

(<!>) = {ALL)-''^^A'h{0,A,MALLY) , (5.3a) 
with the universal scaling function 



r(y/2) VV2 y, 



y 00 
(5.3b) 



For analyzing the numerical data, it is useful to define 
the ratios 



V 



($2 



1, 5"= 1 



3($)($2) 



Q = l 



($4 



3($2)2 • 

(5.4) 

Note that the ratio Q is identical in form to the well 
known Binder cumulant for equilibrium systems. From 



CP 
sDP 

2(1-1/x)/3 



d=5, pbc 



CP 
sDP 

2"'r(.v/2+1/2)/r(A72) 
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Figure 2: The universal order parameter scaling function Mi 
(inset) and the universal fourth order ratio scaling function 
Q as a function of the rescaled source. 

Eqs. (5.1, 5.2), we immediately obtain 



V 



yr(y/2)' 



2T{iy + l)/2y 




y^o, 



(5.5a) 



S = 



A 2vJ 



(5.5b) 



f)- o-K'-'' 

2yJ 3V 

with the scaling argument y = Ahh{ALL)'^. 

Figure 2 compares our analytic results for the normal- 
ized order parameter Mi and the ratio Q to our numerical 
findings. The solid dot marks the condition Q = for 
y = I, and the horizontal dashed line corresponds to the 
limit 2/3. Fig. 2 demonstrates that the data of the lat- 
tice models obey the modified FSS form Eqs. (5.3b) and 
(5.5b), and that the obtained scaling curves are in perfect 
agreement with the results of the continuum theory. We 
rate this as an impressive manifestation of the robust- 
ness of the DP universality class. Two further points 
are worth stressing: (i) The order parameter assumes 
both asymptotic regimes (Mi ~ ^/y for y —> 00 and 
Ml ~ tt/2 y for y ^ 0) predicted by our theory, (ii) 
As mentioned above, the simulated systems got stuck 
quickly in the absorbing state if the external source was 
turned off, h = 0. Thus both, the analytical results as 
well as the numerical simulations refiect that well-defined 
steady-states exist close to the critical point for > 
only. 

The corrections due to the higher modes {v > 0) and 
the exponential instead of the algebraic crossover to mean 
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field scaling for y ^ oo are not resolved by the numer- 
ical data. Note, however, that the leading terms of the 
order parameter and second moment, cf. Eqs. (5.3b) and 
(5.5a) are correct in this limit. In this mean field region 
the order parameter fluctuations are dominated by small 
Gaussian correlations. Hence, we have ($^)/(4>)^ = 
1 + {N{N - 1)/2)V + 0{V'^). Using this expansion one 
easily demonstrates that S and Q as given in Eq. (5.5b) 
show the correct asymptotic scaling including the correc- 
tions cx y~^- 

In the absorbing state, the ratios Eq. (5.5) are not finite 
for y 0. To analyze the scaling behavior in this limit, 
we introduce the following combination of moments: 



U = 



2-3S 
2-3Q 



(5.6) 



which can be viewed as an analog in critical dynamics of 
the famous Binder cumulant. Inserting the lowest-mode 
scaling functions, Eq. (5.5b), this ratio becomes simply a 
constant equal to 1/2. This value is indeed correct in the 
limit y ^ oo, but for y ^ 0, we should expect deviations 
due to the finite-size shift of the critical control parameter 
T. Using the scaling form (4.18) for the order parameter 
moments, we obtain 



C/(a,6) 



di{a,b)di{a,b) ~ di{a,b)d2{a,by 



(5.7) 



as functions of the parameters a and b, with i?;(a, b) fol- 
lowing from Eqs. (4.15) and (4.17b) via the recursion re- 
lation (4.16). U{a,b) is displayed in Fig. 3 as a function 
of In 6 with a as a parameter. It is only in the case a = 
that U is constantly equal to 1/2, whereas U deviates and 
grows from this value for 6^1 when a gets increasingly 
negative. Thus, the limit Uo{a) = lim&^o C^(a, ^) at the 
bulk critical point is a convenient measure of the shift of 
the critical control parameter due to finite size. 

Now, let us look more closely at the parameter a as 
given in Eq. (4.7a). For x = y = wc get 



l_Z!!!7^(2)(o) 
4 



0(w*)] {l + O(i;2)}vi?(i)(0). 

(5.8) 

Defining a new nonuniversal length Lq by the relation 



vD^'\0) = iLo/Lf/' 



(5.9) 



we obtain 



1 - 



47:>(i)(0)2 
0{{Lo/V''-' 



(Lo/Lf-^ + OiiLo/Lf"-') 
)}iLo/Lf/^-\ (5.10) 




log b 



Figure 3: The ratio U as a function of In 6 with a as a pa- 
rameter. For a — 0, U lies on the abscissa, U = 1/2. For 
a decreasing form zero, the values of U deviate increasingly 
from 1/2 for 6 < 1. 



At first glance, the correction factor 
-715(2) (o)/4i:)(i) (0)2 ^ -2.096 in d = 5 seems to 
be a universal contribution. Note, however, that this 
correction factor merely represents the 1-loop contri- 
bution and that the 2-loop contribution 0{{Lo/L)'^~'^) 
of the second factor in Eq. (5.10) is of the same order 
in Lq/L as the 1-loop contribution. Therefore, to be 
consistent, one cither has to take only the lowest order in 
Eq. (5.10), or, if one seeks to proceed to next to leading 
order, one has to account for the 2-loop contribution 
to the finite-size shift of the control parameter r. This 
subtlety was overlooked by Chen and Dohm [10] in their 
work on EES in the Ising model, and their derivation of 
universal scaling functions. Because there exists to date 
no 2-loop calculation of the shift of the critical control 
parameter, which would eventually lead to a universal 
correction proportional to {Lq/ L)'^^'^, we introduce 
here a wild-card K for this universal correction. The 
introduction of K produces 

a = -[l + K {Lo/Lf-^ + OiiLo/L)^'-')] {Lo/Lf/^-' , 

(5.11) 

with Lq and K to be determined by fits to the numerical 
data. 

Next, we revisit the ratio U. From the representa- 
tions (4.16) (4.17b) of the functions i?; in Eq. (5.7), we 
obtain in the limit b ^ 



U = f7o(a) := U{a,0) 

_ [F{a)~a] [{l + a^)-aF{a)] 
~ (2 + a^)F{a) - a{3 + a^) 



(5.12a) 
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Figure 4: (Color online) The ratio U as a function of L. The 
meaning of the various curves is explained in the text. 

where 

F{a) = y27^cxp(-aV2) crfc(a/V2)~\ (5.12b) 

with a = y^Lo/L (l + KLq/L + 0{{La/Lf)) . The ratio 
U is shown for d = 5 in Fig. 4. The solid dots stem 
from our Monte-Carlo simulations of critical sDP on bcc 
lattices of linear size i = 4 to i = 32. The red up- 
per curve is a fit to the numerical data with Lq = 1.01 
and K — 2.17. As expected, the nonuniversal length 
scale Lq is of the order of the lattice spacing. Note that 
the numerical result for correction parameter K is posi- 
tive whereas the pure 1-loop calculation, which entirely 
misses the 0((Lo/i)''~'*)-tcrm in Eq. (5.10), pretends a 
negative value of —2.096. For the purpose of demonstra- 
tion, we include in Fig. 4 the curve of Uq pertaining to 
this fiawed value of green lower curve. Note that this 
curve has a non-physical maximum near L = 6, which 
clearly shows that a pure 1-loop calculation is incomplete 
and misleading and which underscores our previous rea- 
soning that 2-loop contributions to the shift of the control 
parameter cannot be neglected for the interpretation of 
the numerical data. For further comparison, we also plot 
[/q using the correct 1-loop result for the control parame- 
ter shift, i.e, with K set to zero (blue middle curve). The 
figure shows that the corrections cannot be neglected be- 
low L w 24 due to the slow decrease of a ~ ^L^jL with 



increasing L. The zero-mode limit 1/2 (brown straight 
line) approximates V reliably for only very large values 
of L, which were out of reach for our simulations. 

VI. CROSSOVER TO MEAN FIELD BEHAVIOR 

In this section, we consider the crossover to the mean 
field behavior in the infinite volume limit w — s- cx3 in spa- 
tial dimensions d > 4. If w is comparable with or greater 
than 1, i.e., for ?■ = (t — Tc) -I- M > {2tt/LY , we can and 
do calculate the order parameter M = (<I>) = g{s{r,t)) 
and its correlation x'^ ~ {S^'^) ~ (<I>^) — (<I>)^ = 
g'^L-'^J d'*r(s(r,t)s(0,t))=™^ in a standard 1-loop per- 
turbation expansion based on functional integrals with 
weight exp(— >/), where is the response functional as 
given in Eq. (2.3). The results for M and produced 
by this direct calculation are then compared with the cor- 
responding expressions calculated with the steady state 
distribution function, Eq. (4.12). This comparison re- 
veals that neither the lowest-mode approximation nor 
the 1-loop calculation using the Markovian approxima- 
tion capture the correct crossover behavior as produced 
by the direct calculation. Equation (6.9) nails down the 
difference in the outcome of the direct calculation and 
the one that uses the Markovian approximation. 

For calculating M and without recourse to the 
dynamic free energy, we need to know both the propaga- 
tor and the correlator implied in J7. Applying the shift 
s ^ s + M / g to j/, we readily obtain 

G(q, t) = e{t) exp(-A(r + q2)t) , (6.1a) 

g(q,0 = ^(^^^q2) g(q^l^l)^ (6.1b) 

as the propagator and correlator in time and momentum 
space. Then, to 1-loop order, the equation of state follows 
as 

= rAf+iAf2 + ^L-rf^C(q,0). (6.2) 

q 

To the same order we obtain for the correlation 
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^X"' = C{Q, 0) + {XgfL-" T.JJ'^^ [^(0' ~t)G{0, (^C^lq, t - t' f - 2C(q, t - t')G(q, t - t') 

^ — oo 

+ G(0, -t)C(0, -t') (2C(q, t - t')G(q, i - i') - G(q, t - t')^ 
I 



(6.3) 



After some calculation and after rescaling of the param- 
eters and fields as before, we get the equation of state, 



y = {x + z)(\ + 



l-wD'^^\w)\z z2 



In the active region —x ^1,2;^::^ 4j/, we get 
M ~\x\[l + ^ 



(6.4) 



in terms of the scaling variables defined by Eqs. (4.1) 
and (4.2). The correlations in terms of these variables 
are given by 
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—X ' ~ 1 - — 
g-^ x'^ 



(6.7a) 



(6.7b) 



2L'' _i _ z 
g'^ X + z 



2 - 4u;D(i)(w) + 3w£)(i)(3u;/2) Finally, we find for 4?/ > 



+ z 



{x + zY 

[x + z)3 



(6.5) 



If w 3> 1, we have D'^^\w) 



up to exponentially 



small corrections. Using these properties, we find that 
the equation of state and the correlation approach their 
mean field forms in the infinite volume limit with expo- 
nentially decreasing deviations. In contrast to this expo- 
nential crossover, the lowest-mode approximation, which 
corresponds to letting D^^\w) — > 0, produces unphys- 
ical algebraic crossover to mean field behavior with de- 
creasing deviations proportional to (a; -I- z)~'^. 

Recall that we have calculated in Sec. IV with the help 
of the steady state distribution a scaling form for the 
moments of $, Eq. (4.18). This equation implies scaling 
forms for the equation of state and the correlations, which 
we in the following wish to compare to Eqs. (6.4) and 
(6.5). For simplicity, we focus on the following three 
regions of phase space: the absorbing phase region x ':$> \, 
the active phase region ~x 3> 1, both with small source 
x"^ 2> 4y, and the region with large source including the 
bulk critical point 4y ^ x^. Expanding Eqs. (6.4) and 
(6.5) for a; ^ 1, x^ ^ Ay, we obtain 



id/2 



-M 



1 



x^ 



r 



x'^ \ x"^ 



3 - 5wL»(i)(«7) -I- 3u;i:'(i)(3u;/2) 



(6.6a) 



(6.6b) 



M~V^ 1- — z-^h (6-8a) 



9 
1/ 



2^/y 



4y 



1 X 



-X 



2 4Vy 



1 - AwD'-^^ (w) + SwD'^'^^ {3w/2) + w'^D^'^^ {w) 

(6.8b) 

Next, let us see what our steady state distribution, 
Eq. (4.12), tells us, and let us compare that to the above. 
Using the asymptotic properties of the parabolic cylin- 
der functions [36] in the three regions, we obtain for M 
the same expressions as displayed in Eqs. (6.6a), (6.7a) 
and (6.8a). For the correlations, we recover Eqs. (6.6b), 
(6.7b) and (6.8b) up to one alteration: the function 
D^^\?>'w/2) is replaced by 

D^^\-iw/2) ^ D'^^\w) - ^D^^Hw) , (6.9) 



which is an identity to linear order in w but which mod- 
ifies the correlations at higher orders. The mean-field 
parts of the expressions for given by the respective 
first two terms on the right hand sides of Eqs. (6.6b), 
(6.7b) and (6.8b), are identical in both approaches. For 
w ^ 1, where we have D^^^ (w) ~ w^^ up to exponentially 
small corrections, Eqs. (6.6) to (6.8) tend to the mean- 
field behavior with exponentially decaying corrections. 
After the replacement (6.9) (i.e., in the approach based 
on the steady state distribution), however, these correc- 
tions for the correlation fall off only algebraically. 
This incorrect feature is a consequence of the Markovian 
approximation as the direct calculation shows. 
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VII. FINITE SIZE EFFECTS BELOW dc 

As mentioned above, in a former publication [17], one 
of us and coworkers have calculated finite size effects 
for absorbing nonequilibrium processes belonging the DP 
universality class in spatial dimensions cZ = 4 — e < 4. 
There, systems without a source were considered, and the 
consequences of the finite size scaling for the relaxation 
behavior were scrutinized. Here, we are interested in the 
steady state properties in the presence of the source h. 
We calculate various quantities in an e-expansion, most 
notably the parameters a and b and the ratio U. Equa- 
tion (7.17) gives our e-expansion results for a and 6, and 
Eq. (7.18) states our e-expansion result for U. 

A 1-loop calculation for d < 4 can be done in much 
the same way as the calculation for d > 4 presented 
in Sec. III. Now, however, in addition to the functions 
£)(')(u!) in Eqs. (3.14), the brackets in Eqs. (3.13) become 
IR divergent, and therefore, they no longer can be simply 
included in non-universal amplitudes. Rather, these ad- 
ditional divergencies must be handled with a renornial- 
ization procedure and the renormalization group equa- 
tion. For general background on these methods, we re- 
fer to the usual textbooks, e.g. [8, 41]; for applications 
of these techniques to the DP universality class consult, 
e.g., Refs. [16, 23, 24]. 

Whereas the coupling constant g is dangerously irrel- 
evant in d > 4, it is, respectively, marginal and relevant 
in d = 4 and d < 4. Therefore, it is useful for the case 
d < 4 presented in this section and the case d = 4 to 
be presented in Sec. VIII to recast the dynamic free en- 
ergy (3.5) as 



\-^kdt + f + ^{S~S) S-hSy 

(7.1) 

with S{t) = L^"^ J d'^r sir, t) and likewise for s. To facil- 
itate the renormalization procedure and to cleanly keep 
track of bare (unrenormalized) and renormalized quanti- 
ties, we henceforth label bare fields and parameters with 
a ring °, i.e., we let s — > s, s ^ s, t — > f, and so on, and 
we reserve symbols without a ring ° for their renormalized 
counterparts. The bare and the renormalized quantities 
are related via the renormalization scheme 



T = ZZ^ZrT - 



s = Z^''^s, h^Z-^Z^/^h, (7.2a) 
g^ = ZfZ-^Z^g^ , (7.2b) 
Tc . (7.2c) 



mined as to eliminate the e-poles arising in a dimen- 
sional regularized calculation of the momentum space 
integrals. This kind of calculation orders naturally in 
powers of a dimensionless coupling constant u defined by 
u = GeH^^g^, where is a convenient length scale, 
and Ge = r(l + £/2)/(47r)''/2. The renormalization fac- 
tors are given to 1-loop by 



Z =1 

Zr^l 



u 

u 
2e 



Zx = 1 
Zu = 1 



8e ' 
2ii 



(7.3a) 
(7.3b) 



With help of the renormalization scheme (7.2) and the 
renormalization factors (7.3), we find that the renormal- 
ized versions of the parameter functions featured in the 
dynamic free energy (7.1) are given by 



2 V 27r / 4 ^ . 

u 



27r\2 



L 



<7{w) — wa'{w) 



Here we have defined the function 

a{w) = w{\nw- 1) - \d^^\w) , 



(7.4a) 

(7.4b) 
(7.4c) 

(7.5) 



The renormalization factors Z, Z-r and so on are deter- 



where it is understood that _D'^' {w) is taken at d = 4 and 
where a'{w) stands, as usual, for the derivative of <7{w). 
The virtue of the function a{w) is that it and its deriva- 
tive lack the non-analytic behavior of D'^^^w) for w — s- 0. 
However, as shown in the appendix, these functions are 
nevertheless logarithmically divergent in the bulk limit 
w — > oo. In principle, one should handle these diver- 
gences by subtracting a term wln(l -|- w) or ln(l -I- w), 
respectively, as done in [17]. These subtractions, with w 
as given by Eq. (3.11), combine with the logarithm in the 
first brackets of Eqs. (7.4) to produce the IR-divergent 
term \n{{2TT / fiL)^ + {t + M)/fj.'^), which should be elim- 
inated by the renormalization fiow. Nonetheless, we can 
here set these subtleties aside and ignore the divergences 
for w ~> oo because we arc only interested in the strong 
finite size case w 1. 

The perturbation results for the parameter functions, 
Eqs. (7.4), cannot be used directly as they stand. These 
results must be transported by the renormalization group 
fiow to a non-critical region. To this end, we derive in a 
standard fashion Gell-Mann-Low renormalization group 
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equations (RGEs) for the parameter functions via ex- 
ploiting the fact that the bare theory must be indepen- 
dent of the length scale introduced by renormaliza- 
tion: 



In fc = 7 , 
Pplnf = 7 - C, 
2?^ln.9 = i(37-2C). 



(7.6a) 
(7.6b) 
(7.6c) 



where, stands for the renormalization group differen- 
tial operator 

= ^J.^fJ, + (3du + AC^A + TKdr + "y ~ u ~ ^) ' 

(7.7) 

and where 7, Q and so on are the usual RG functions. For 
DP, these RG functions are known to 2-loop order [23, 
24]: 



7 = -- + f2-31n. , 
' 4 V 3/32 



c = - 



8 



17- 2 In- 



4\ 



3/ 256 



'^=¥-(^+10^^^3)256' 

3ii^ / 4\ 
/3 = -s.+ — - (169 + 106 In-) — 



(7.8a) 
(7.8b) 
(7.8c) 
(7.8d) 



where we have included the 2-loop contributions, even 
though we work in this section only to 1-loop order, be- 
cause we will need them in Sec. VIII. The RGEs can 
be solved using the method of characteristics. The idea 
behind this method is to consider all the scaling param- 
eters as a function of a single flow parameter i. One sets 
up characteristic equations that describe how the scaling 
parameters transform under a change of L The charac- 
teristic for the inverse length scale fi is trivial and has 
the solution /2(£) = i.e., a change of I corresponds to 
a change of the external inverse length scale. With help 
of the solution to the remaining characteristics and also 
with help of a dimensional analysis to account for naive 
dimensions, we obtain 



k{T,M,u,fi,L) = Xiey^k 



--{£) M{e) 



,u{e),l,fj.£L 



(7.9a) 



„ / f{e) M{e) 



f{e) M{£) 



,u(^),l,Mtt , (7.9b) 



X g 



,u{e),i,fieL , (7.9c) 



where 

f{e)=TXr{i), (7.10a) 
M{£) = M£'^^[u{e)/u]^/^X{£)-^/^. (7.10b) 

At this stage, the scaling relations (7.9) are still rather 
formal because we still must determine X{£), X\{£), 
Xr{£) and u{£) by solving their respective characteris- 
tic. The characteristic for the dimensionless coupling 
constant u is given by 



(7.11) 



where we abbreviated v = u{£). The remaining charac- 
teristics are all of the same structure: 



' d£ 



q{v). 



(7.12) 



Here, Q is a placeholder for X, X^, and X\, respectively, 
and g is a placeholder for 7, k, and Q, respectively. As 
usual, solving the characteristics leads to qualitatively 
different results depending on whether we consider the 
upper critical dimension or dimensions below it. We will 
return to d = 4 in Sec. VIII. 

For d < 4, the dimensionless coupling constant u flows 
to the stable fixed point u^, = 2£/3 -I- O(e^), and, conse- 
quently, X{£) etc. display power law behavior described 
by the well known critical exponents of the DP univer- 
sality class. Using a compact notation where p stands 
ambiguously for k, f and we can write the resulting 
scaling form for the parameter functions as 



p(t, Af, L) = £^fp{£-^^''T, £-l^'''M, £L) , 



(7.13) 



with (5^ 



-77 = d — 2/3/i^, 5f = z — T] = ^jv and 6g = 
{2z~d — 3rj)/2 ~ (7 — /?)/;/, respectively. Here, the three 
independent critical exponents are given by 



6 



7 = l+g, 



1 e . X 

1^=- + — , 7.14) 

2 16 ' ^ ^ 



up to terms of order [23, 24]. The exponents /3 
and 7 must not be confused with the RG functions 
discussed above. Now, we choose the flow parameter 
£ = 2tt/i^iL <C 1 to eliminate the IR-diverging logarithm 
ln(/^L/27r). The parameter £ must be small to reach the 
asymptotic region, i.e., to produce universal behavior. Of 
course, this is a condition on the size L, which must not 
be small in comparison to a non-universal length scale 
Lq which is set in our simulations by the lattice con- 
stant. After implementing our choice of £, we obtain the 
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basic parameter functions in scaling form 

.2f)/v-d 



0.9 



k = {fj,L/2TT 



1 + -a'{w) 



l + -a'{w)\T{^lL/2^) 
e 
6 



a{w) ~ W(j' {w) /i^l 



9 = (/.iL/27r 
where now 



,(/3-7)/!' 



l + -<j'{w) 



(7.15a) 

(7.15b) 
(7.15c) 



[T{fiL/2TTY^" + M{nL/2TTf'"]/fi^ . (7.16) 

Next, let us return to the parameters a and b of the 
dynamic free energy (4.5) and the steady state distribu- 
tion (4.12). Because we are interested in the strong fi- 
nite size case if ^ 1, we can approximate a{w) w cr(0) + 
w'(0) and a'{w) « cr'(O), where cr(0) = -81n2/7r2 « 
-0.56184 and cr'(O) = -I - Ce - 2 ln2/3 - 6C'(2)/7r2 « 
— 1.85789 with Ce and C denoting Eulcr's constant and 
Riemann's C-function, respectively. Recalling the defini- 
tions of a and b, Eqs. (4.6), we find after some algebra 
their e-cxpansions to be given by 

a =7rv/67^ [1 - ^£ + 0{e^)] {t / i?){piL /2ttY''' 

- ^yi76[l + 0(£)] , (7.17a) 

TT 

b ^ (SttVe) [1 - 2Ae + 0(e2)] (i///)(^L/2^)^/'^ , 

(7.17b) 



where A = (ln7r)/4 - (ln2)/12 + [Ce - 1)/8 - 
3C'(2)/(27r)2 « 0.24688 and A = /3 + 7. Now, we are 
finally in the position to address our main observable, 
the momenta ratio J7o(a) = U{a,b 0). Expanding 
Eqs. (5.12) with a as given in Eq. (7.17) in e wc find that 



Stt Vtt 



2e /, 4 16 \ , 



0(e3/2) 



(7.18) 



at the bulk critical point t = 0. Setting e = 1. 2, and 3 
in Eq. (7.18) we obtain estimates of Uq for systems be- 
longing to the DP universality class in spatial dimensions 
3, 2, and 1, respectively: 



Un 



These results are to be compared with our numerical 
data for d below dc ~ 4. For d < 4, we simulated in 
addition to the CP and sDP also the pair contact pro- 
cess (PCP) introduced by Jensen (see [42] as well as 
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Figure 5: The universal ratio U at the bulk critical point as a 
function of the scaled source h/L^^" in dimensions d = 1, 2, 3. 



the reviews [13] and [15]). The corresponding numer- 
ical curves for U are shown in Fig. (5). From these 
curves we find U{d = 3) = 0.61, U{d 2) 0.704, 
and U{d = 1) = 0.833. Our 1-loop calculation repro- 
duces qualitatively correct the right trend of C/ as a func- 
tion of d. Not surprisingly, the quantitative agreement 
is rather poor for low dimensions. At least for d = 3, 
the e-expansion estimate is not too far away from the 
numerical value. For a field theoretic 1-loop calculation 
of amplitude ratios, errors of about 10% are typical for 
£ = \ [7, 8], and the deviation of our analytical and nu- 
merical results for £ = 1 is consistent with that. It is im- 
portant to note that the value of U (d) for a given dimen- 
sion is quantitatively the same for the three processes, 
and that, therefore. U{d) proves to be a true universal 
signature of the DP class. Moreover, the universality of 
U (d) for the three processes shows that PCP definitely 
belongs to the DP universality class. 



VIII. FINITE SIZE EFFECTS AT dc 

Here, we study finite size effects right at dc, where 
the finite system size is expected to generate logarith- 
mic corrections to the bulk behavior. Guided by lessons 
learned form previous studies of logarithmic correc- 
tions [29, 40, 43, 44], we choose to derive scaling forms in 
a parametric representation rather than in the more tra- 
ditional representation featuring nested logarithms. Cen- 
tral to the parametric representation is the parametriza- 
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tion (8.4) of the system size L. Equation (8.8) gives our 
general parametric results for a and &, and Eq. (8.10) 
specializes the result for a to the critical point. Finally, 
we compare our parametric result for U to our simulation 
results. 

Past studies of logarithmic corrections in DP [29, 40] 
and other systems, e.g., linear polymers [43, 44], led to 
the observation that one has to push the analytic calcula- 
tions beyond the leading logarithmic correction to obtain 
good agreement between theory and simulations. To go 
beyond the leading logarithmic correction, we will work 
in the following, as announced above, to 2-loop order as 
far as the RGEs are concerned. Concerning the scaling 
functions, it will still be sufficient, for the most part, to 
work to 1-loop order. However, here is an important ex- 
ception: a 1-loop calculation of the scaling function of 
a does not suffice to determine the next to leading log- 
arithmic correction to a entirely. This subtlety will be 
discussed as we move along. 

Our vantage point for this section will be the gen- 
eral scaling forms for the parameter functions fc, f and 
g derived in Sec. VII, Eqs. (7.9) in conjunction with 
Eqs. (7.10). To fill these general scaling forms with live 
for d = 4, wc must solve the characteristics for this di- 
mension. In order to make our notation as compact as 
possible, we will write in the following the RG functions 
as f{u) — fo + fiu + f2u'^ + ■ ■ ■ with / standing ambigu- 
ously for 7, C, K, and f3. The meaning of the coefficients 
/o, fi and should be evident. 

First, we solve the characteristics for d ~ 4. The solu- 
tion to the characteristic for the dimensionless coupling 
constant u, differential equation (7.11), is given by 



where Lq = 2tt/{^£o). Note from this relation that, in 
contrast to the 1-loop approximation of the RG functions 
(/^a — ^ 0), the 2-loop approximation leads to an effective 
i-depcndence of the nonuniversal length Lq which must 
not be neglected (see also our discussion of Fig. 6 below). 
Taken together, Eqs. (8.2) and (8.4) can be exploited as 
a parametric representation of the tuple (L, Q) with v as 
parameter. This representation has the advantage that 
the resulting formulas are comparatively compact and, 
more importantly, that one deals with clean expansion 
in powers of v. 

After this prelude, let us return to the parameter func- 
tions. Collecting from the renormalized perturbation cal- 
culation results (7.4), the general scaling forms (7.9), the 
solutions of the characteristics (8.2), and implementing 
our choice (8.3) of the flow parameter, we obtain 



,-1 3 , 



a (w) 



-1/6 



cxp[cf^v + 0{v^)] , (8.5a) 



V ^ a'iw) 
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1 /27r 
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5 = 47r 



[^(y^) _ ^a'{w)] cxp cf'v + 0{v^) 
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exp [cgv + 0{v'^)] , (8.5c) 
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e = £{v) =£ou 



exp 



i/32v)-' + 0{v) , (8.1) 



with £() being an integration constant. The characteris- 
tic (7.12) is readily solved with the result 



(g2/32 - gi/33) 
/3| 



(8.2) 



with a non universal integration constant Qq. 

Next, we choose the flow parameter £ such that the lat- 
tice size L effectively acquires a finite value in the scaling 
limit: 



'^ = 1 
27r 



(8.3) 



With this choice, £ and v tend to zero for /zL c», and 
L and v are related via 



(L/Lo) = v^'/^" exp \{f32v)-' + 0(^)1 , 



(8.4) 



where now 



w 



,1/4 



2tt 
I? 



exp [caV + 0(v^ 



+ m — u^^/^^ cxp 



(8.6) 



and where we have included nonuniversal integration con- 
stants stemming from characteristics solutions (8.2), viz. 
^0- X\,o Siiid Xr.o, in the nonuniversal amplitudes of r, 
h, and to. The coefficients appearing in the exponentials 
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in Eqs. (8.5) and (8.6) are given by 
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(8.7b) 
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(8.7d) 

0.07971 . 
(8.7e) 



Now, we revisit a and b. Inserting our results (8.5) into 
definitions (4.6) we find 



hi 



1/4 



exp [caV + 0{ 



+ '- [v-^ + F{w)\ [u{w) - wa'iw)] exp [0{v^)] 

(8.8a) 
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where 



U ^ cr'(w) 
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exp [c,,u + 0(u2)] , (8.8b) 
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-0.00724. 



(8.9) 



At this point, a comment is in order. A full-fledged 2-loop 
calculation of a's universal scaling function is expected 
to produce, inter alia, terms of the same order in v as 
the 1-loop calculation. Therefore, we have replaced in 
second line of Eq. (8.8a) the 1-loop contribution i;^/^ by 
the bracket containing F(w), where F(w) is a hitherto 
unknown function. We will leave the calculation of F{w), 
which will be challenging, to future work. 

Finally, let us return to our ratio U of the order pa- 
rameter moments. As was the case for d < 4, we are 
mainly interested in the strong finite size regime w 1 
and. therefore, we approximate a(w) « cr(0) -I- wa'(0) 



and (7'{w) w cr'(O). Focussing on criticality, we set r = 0. 
The remains a, Eq. (8.8a), are then 



4 In 2 



K 



-1/2 



(8.10) 



with a universal correction K = F{0). As discussed 
above, a calculation of K would require to determine 
the scaling function of a to 2-loop order. Because corre- 
sponding results are currently not at our disposal, we 
use K as a fitting parameter. Note that a falls off 
only as a [ln(L/Lo)]~^/^ in d = 4 compared to the 
a ~ (L/Lq)"^/^ behavior in c? = 5. Thus, it must be 
expected that U approaches its zero-mode limit 1/2 even 
slower for increasing system size than in d = 5, and that 
one needs at dc even larger systems than above c?c for the 
zero-mode theory to provide a good approximation. 

Substituting Eq. (8.10) without a further expansion 
into our scaling function C/o(a), Eq. (5.12), we obtain our 
final result for Uq as a function of v. We then use the so- 
obtained expression for U in conjunction with Eq. (8.4) 
as a parametric representation of the tuple (L, U) with v 
as parameter, which we plot together with our numerical 
data in Fig. 6. In the plot, we use K and Lq as fitting 
parameters. Our best-fit analytical curve (the red mid- 
dle curve) impressively tracks our data points over the 
entire range of simulated lattice sizes including sizes as 
small as L = 4. As mentioned earlier, the 2-loop RG 
contribution to Eq. (8.4) effectively modifies Lq. Due 
to this modification the slope of the continuous curves 
is significantly reduced for L < 40 in comparison to the 
dashed pure 1-loop curve. The introduction of K leads 
mainly to a rescaling of the nonuniversal length Lq which 
manifests itself in the modest deviation of the red middle 
curve from the blue lower curve for L < 10. Note that 
up to 2-loop order, one can eliminate K entirely from U 
via a simple rescaling of Lq. Thus, one may view the 
introduction of a non-zero K as a crude way of account- 
ing for the infiuence of high loop-orders. Note also, that 
the effect of the non-zero, fitted K is much smaller than 
the effective modification of Lq resulting from the 2-loop 
RG contribution to Eq. (8.4). Over all, the agreement 
between theory and simulation is remarkable. This ob- 
servation reassures us once more about the validity of 
our analytical and numerical approaches. Moreover, it 
underscores the advantages of the parametric representa- 
tion and makes tangible the necessity of including 2-loop 
RG results. 
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Figure 6: (Color online) The universal ratio U versus system 
size L for d = 4. The blue lower and the red middle curves 
represent our analytical results (5.12) and (8.10) with Lq — 
2.4, K = and Lq = 1.5, K = —0.5, respectively. For 
comparison, we included the dashed upper curve, where we 
have disregarded any 2-loop contributions and where we have 
fitted Lo to the data points for larger L, La = 5.6. The solid 
dots stem from our Monte-Carlo simulations of critical sDP 
on bcc lattices of linear sizes ranging from L = 4 to L = 64. 

IX. CONCLUDING REMARKS 

In summary, we have investigated finite size scaling 
effects in steady state systems belonging to the directed 
percolation universaUty class. We have assumed a hyper- 
cubic geometry with length L, periodic boundary con- 
ditions and the presence of an external homogeneous 
time independent source which prevents the systems to 
fall into their absorbing inactive state. We applied a 
field-theoretic technique based on an effective response 
functional (dynamic free energy) for the lowest (homoge- 
neous) mode, which allowed us to calculate finite size ef- 
fects within a 1-loop perturbation expansion of the higher 
modes combined with a Markovian approximation. This 
latter approximation is indispensable for calculations of 
strict non-equilibrium properties of systems without de- 
tailed balance. In particular, it allowed us to calculate 
the steady state distribution for the lowest mode via the 
associated Fokker-Planck equation. Using this distribu- 
tion, we calculated explicit scaling forms for the moments 
of the homogeneous order parameter. Moreover, we in- 
troduced and calculated a ratio U of order parameter 
moments which allowed us to analyze universal finite size 
effects right at the critical point. Complementary to our 
analytical work, we performed Monte Carlo simulations 
based on the contact process, the site directed percola- 
tion process and, on occasion, the pair contact process. 



Above and at the upper critical dimension 4, we found 
remarkable agreement between our analytical and numer- 
ical approaches. In these dimensions, the usual coupling 
constant of the cubic term in the response functional 
is dangerously irrelevant. Due to this dangerous irrel- 
evance, the universal scaling functions depend on the 
additional (compared to c? < 4) scaling variable L/Lq, 
where Lq is a nonuniversal length scale. Our results 
demonstrate that it is necessary to push the diagram- 
matic calculations beyond 1-loop order to obtain agree- 
ment between theory and simulations down to very small 
systems sizes, L/Lq ~ 1. 

For d below 4, we calculated the universal critical val- 
ues of ?7 in a e-expansion to order 0{e^^'^). The accu- 
racy of this calculation corresponds to that of the calcu- 
lation of the Binder cumulant of the c/)*-model at the bulk 
critical point by Brezin and Zinn- Justin [7]. The agree- 
ment between our theory and simulations is within the 
expectation for a 1-loop calculation that captures terms 
to 0(£'^/^), and, of course, it decreases for decreasing 
dimensions. However, the universal critical values of U 
produced by our simulations were identical for all three 
processes that we simulated and, therefore, U proved to 
be a true signature of the DP universality class. More- 
over, this finding demonstrates that the pair contact pro- 
cess belongs to this class. 

On the analytical side, our study is the first investiga- 
tion, besides the former work of one of us and coworkers, 
which addresses finite size scaling near absorbing phase 
transitions. We believe that our approach may be applied 
to many other non-equilibrium phenomena, and that it 
can help to improve the understanding of finite-size ef- 
fects in non-equilibrium systems significantly. 

Appendix: PROPERTIES OF THE FUNCTIONS 

Ll('^(w) AND a{w) 

Using an exponential representation of the de- 
nominators in the sums (3.7), (r -I- q^)~' = 
r(Z)~^ /o°° dtt''^^ exp(— (r + q^)t), we eventually obtain 
for the functions 13*^'' (w) the Laplace-transforms stated 
in Eq. (3.12). D'-^^w) and D^^^w) arc smooth func- 
tions if d > 4 (with D(i)(0) = 4.229 and = 
21.421 for d = 5). Due to the recursion relation 
£)('+!) (■^^ — QD^'-) {w) / dw, we can restrict our atten- 
tion here to £'(^'(w) in order to determine the remaining 
properties of D*-''(i(;) that are used in the main text. 

To extract the behavior of D^^^ [w) at small arguments. 
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we divide it into parts, 

D^^\w) = Ii{w)+l2{w) 

with 
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where k is some integer with k > d/2 — 1 to provide 
intcgrabihty at t = 0. The specifics of the remaining part 
l2{w) can easily be gathered from Eqs. (3.12) and (A.2). 
Then, it is straightforward to see l2{w) is an analytic 
fimction of w. Hence, we obtain for small w = that 



r(l-d/2) 

(A.3) 

where D^^-' [w) is analytic. 

To extract the behavior of Z?^"'^-' {w) for large arguments, 
we divide this function in three parts 

D^^\w) = .h{w) ^ J2{w) + .h{w) . (A.4) 

The behavior of Ji(w) and J2{w) for w 3> 1 is given by 

d/2 



Ji{w) = 

J2{w) = 



dte-""* 



(f) -A(t)'^ + l] =0(c-™), 

(A.5) 



Using the expansion A{t)'^ — 1 
0(exp(-2t)), we find 



(A.6) 

4dexp(-t) + 



dte 
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- A{tf 



dsC-^''n,/s 



d/2-2 r 



1 -^(s)-^ 

4d / ds { — ] cxp ( — s — Tx^w I s ^ 
%ixdw''l^-^I^Kdl2(2Tx^), (A.7) 



ds ^— J cxp 



where Ka{z) is the Basset function, for the leading be- 
havior of J^lw). Using the asymptotic properties of this 
function, we finally get 

D(i)(u;) = i_47rd«;('^-3)/4 exp(~27rV^) + . . . , (A.8) 
w 

where the ellipsis denote subleading terms. 

For d < 4, we use instead of D^^\w) the function 

a{w) = w{\n w-1)- \d^'^^ {w) , (A.9) 



where _D'"'^-'(w) is taken at d = 4, and its first derivative 
a'{w) to eliminate the nonanalytic logarithmic behavior 
near u; = 0. This function has a power expansion in w 
as derived in [17]: 



(A.IO) 



k=0 



where 



;ln2 



n r ,1 , 21n2 6C'(2) 
CTi = C_E + 1 H 1 ^ 



(A.lla) 
(A.Ub) 



and for k > 2 



(1-1/4'=) 



C(fc)C(fc + i), 



(A.llc) 



with Ce « 0.577716 and C'(2) « -0.937548. For w -> 
oo, the function a{w) behaves as 

1 16 

a(w) ~ w(lnu' — 1) ^ 1 w^^'^ exp(— 27rv^) + . . . , 

(A.12) 

up to subleading terms. 
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